official_birth_records = official_birth_records %>% mutate(births_thousands = births/1000)
official_births_summary = official_birth_records %>%
group_by(country_area, country_area_col) %>%
summarize(min_births = min(births_thousands),max_births = max(births_thousands),mean_births = mean(births_thousands),
max_date = max(date))
ref_date = as.Date("2020-01-01")
g_births = ggplot(official_birth_records, aes(x = date, y = births_thousands, col = country_area_col))
g_births = g_births +
#### reference for variations
geom_segment(data = official_births_summary, aes(x = ref_date, xend = ref_date, y = mean_births, yend = mean_births-5*mean_births/100),
size = 1.5, alpha = 0.7)+
# geom_text(data = official_births_summary, aes(x = ref_date, y = min_births-5*mean_births/100), label = "5% change from the mean",
# hjust = 0, vjust = 0, nudge_x = 100)+
#### uncorrected births
geom_line(aes(y = births_original_numbers/1000), col = "gray80")+
#### corrected births
geom_line()+ #geom_point(size = 0.5)+
#### COUNTRIES
geom_text(data = official_births_summary, aes(x = ref_date, y = max_births *1.05, label = country_area, fontface = 2), vjust = 0, hjust = 1)+
#### settings
scale_color_identity()+
scale_x_date(date_minor_breaks = "1 year", date_labels = "%Y")+
scale_y_continuous(expand = expansion(mult = c(0,0.2)))+ # sec.axis = sec_axis(~ (. - max(.))/max(.)*100, breaks = seq(-30,30,by = 10), name = "% change from the max")
ylab("Births (monthly in thousands)")+ xlab("Date")+
facet_grid(country_area ~., scale = "free")+
theme(strip.text.y = element_blank())
g_births### change
users = read_feather(path = paste0(IO$output_clue, "users.feather"))
n_users_per_country = users %>% group_by(country_area) %>% summarise(n_users = n()) %>%
mutate(country_area_col = dict$country_area$country_area_col[match(country_area, dict$country_area$country_area)],
country_area = factor(country_area, levels = dict$country_area$country_area))
g_n_users = ggplot(n_users_per_country, aes(x = country_area, y = n_users/1000, fill = country_area_col))+
coord_flip()+
geom_bar(stat = "identity")+
scale_fill_identity()+
xlab("")+ylab("K# of app users")+
theme(strip.text = element_blank(),
axis.text.y = element_blank())+
facet_grid(country_area ~ ., scale = "free")
g_n_userspp = read_feather(path = "../Data/outputs/Rdata/aggregated_sex_counts_clue_July2017-June2019_incl.feather")
pp_sum = pp %>%
filter(BC == "all", age_cat == "all") %>%
mutate(sex_type_str =
case_when(
sex_type == "all_sex" ~ "Total Sex",
sex_type == "prot_sex" ~ "Protected Sex",
sex_type == "unprot_sex" ~ "Unprotected Sex")) %>%
group_by(date, sex_type, sex_type_str) %>%
summarize(n_users = sum(n_users),
n_any = sum(n_any),
n_log = sum(n),
.groups = "drop") %>%
group_by(sex_type) %>%
mutate(relative_frequency = n_log / n_any,
relative_frequency = relative_frequency/mean(relative_frequency)) %>%
arrange(sex_type , date) %>% ungroup() %>%
mutate(
sex_type_str = sex_type_str %>%
factor(., levels = c("Total Sex","Protected Sex","Unprotected Sex"))
)
g_sex_freq =
ggplot(pp_sum, aes(x = date, y = relative_frequency, col = sex_type_str))+
geom_line()+
ylab("Relative sex frequency (daily)")+xlab("Date")+
#expand_limits(y = 0)+
scale_color_manual(name = "",values = c("gray40","steelblue2","yellow3"))+
facet_grid(sex_type_str ~ .)+
guides(col = FALSE) +
theme(legend.position = "top",
legend.direction = "horizontal",
legend.title = element_blank(),
legend.background = element_rect(fill = "white", color = "transparent", size = 0.5))
g_sex_freq# import png
clue_screen = image_read("../Figures Tables Media/Media/Clue_app_screens.png")
clue_screen_ratio = image_info(clue_screen)$height / image_info(clue_screen)$width
padding = 20
g_clue_screen = ggplot()+ coord_fixed(ratio = clue_screen_ratio)+
background_image(clue_screen)+theme(plot.margin = margin(t=padding, l=padding, r=padding, b=padding, unit = "pt"))fig_1 = arrangeGrob(
g_births,
g_clue_screen,
g_sex_freq,
ncol = 2, nrow = 2,
widths = c(1.5, 1),
heights = c(0.55,0.45),
layout_matrix = cbind(c(1,1), c(3,4)))
fig_1 = as_ggplot(fig_1) +
draw_plot_label(label = letters[1:3], size = 15,
x = c(0, 0.6, 0.6), y = c(1, 1, 0.45))
fig_1## Loading objects:
## sex_models_df
plotlist = list()
n_items = c()
cas = unique(sex_models_df$country_area)
ok = foreach(ca = cas) %do% {
cat(ca, "\n")
# retrieving the model for this country and only plotting for BC all and sex_type all
j = which(
(sex_models_df$country_area == ca) &
(sex_models_df$BC == "all") &
(sex_models_df$age_cat == "all") &
(sex_models_df$sex_type == "unprot_sex"))
load(file = str_c(IO$p_outputs,"Rdata/sex_models_",ca,".Rdata"), verbose = TRUE)
j = which((sapply(sex_models_this_location, "[[","BC") == "all") &
(sapply(sex_models_this_location, "[[","age_cat") == "all") &
(sapply(sex_models_this_location, "[[","sex_type") == "unprot_sex"))
glm_sex_behavior = sex_models_this_location[[j]]$model
# glm_sex_behavior = sex_models[[j]]$model
# Visualisation of the coefficient
g_coef = ggplot_sex_activity_glm_coefficient(
model = glm_sex_behavior,
show_intercept = FALSE,
show_weekly_patterns = FALSE,
ylim = c(-0.4, 1))
g_coef = g_coef +
ggtitle(ca)+
theme(axis.text.x = element_text(size = 7))
print(g_coef)
plotlist[[ca]] = g_coef
n_items = c(n_items, nrow(g_coef$data))
}## Brazil - Central-West
## Loading objects:
## sex_models_this_location
## Brazil - Northeast
## Loading objects:
## sex_models_this_location
## France
## Loading objects:
## sex_models_this_location
## United Kingdom
## Loading objects:
## sex_models_this_location
## United States - California
## Loading objects:
## sex_models_this_location
## United States - Northeast
## Loading objects:
## sex_models_this_location
fig_2 = ggarrange(
ggarrange(
plotlist[[which(cas == "Brazil - Central-West")]],
plotlist[[which(cas == "United States - Northeast")]],
labels = c("a","b"),
ncol = 2, nrow = 1,
widths = n_items[which(cas %in% c("Brazil - Central-West","United States - Northeast"))] /
sum(n_items[which(cas %in% c("Brazil - Central-West","United States - Northeast"))])
),
ggarrange(
plotlist[[which(cas == "France")]],
plotlist[[which(cas == "United Kingdom")]],
labels = c("c","d"),
ncol = 2, nrow = 1,
widths = n_items[which(cas %in% c("France","United Kingdom"))] /
sum(n_items[which(cas %in% c("France","United Kingdom"))])
),
nrow = 2, ncol = 1
)
fig_2# t = seq(0,4*pi, by = pi/20)
# f = sin(t)
# plot(t, f, type = "l")
t = seq(30*7, 45*7)
d = foreach(ca = G_table$country_area, .combine = bind_rows) %do% {
this_ca = G_table %>% filter(country_area == ca)
res = this_ca[rep(1,length(t)),]
res = res %>% mutate(t = t, d = dnorm(t, mean = this_ca$G*7, sd = Gsd0), G = this_ca$G*7)
res
}
br = seq(min(t), max(t), by = 14)
d = d %>%
mutate(country_area_wrapped = str_replace(country_area, " - ","\n"),
country_area_col = dict$country_area$country_area_col[match(country_area, dict$country_area$country_area)])
g_d = ggplot(d, aes(x = t, y = d, fill = country_area_col))+
geom_area(alpha = 0.6)+
geom_vline(aes(xintercept = G, col = country_area_col))+
scale_color_identity()+
scale_fill_identity()+
scale_x_continuous(breaks = br , labels = br/7)+
scale_y_continuous(breaks = NULL)+
xlab("G, gestation duration\n(weeks)")+
ylab(expression(paste("Density, d",tau)))+ #expression(paste("Value is ", sigma,",", R^{2},'=0.6'))
facet_grid(country_area_wrapped ~ ., scale = "free")+
theme(strip.text.y = element_text(angle = 0, hjust = 0))
g_dggplot(d, aes(x = t, y = country_area, alpha = d))+
geom_tile()+
scale_y_discrete(limits = rev(levels(G_table$country_area)))+
xlab("G, gestation duration\n(weeks)")+
ylab("Density")+
guides(alpha = FALSE)+
scale_alpha_continuous(range = c(0,1))+
scale_x_continuous(breaks = br , labels = br/7)model = image_read("../Figures Tables Media/Media/Figure 3-01.png")
model_ratio = image_info(model)$height / image_info(model)$width
g_model = ggplot()+ coord_fixed(ratio = model_ratio)+
background_image(model)
models = image_read("../Figures Tables Media/Media/Figure 3-02.png")
models_ratio = image_info(models)$height / image_info(models)$width
g_models = ggplot()+ coord_fixed(ratio = models_ratio)+
background_image(models)fig_3 = ggarrange(
ggarrange(
g_model,
g_d,
labels = c("a","b"),
ncol = 2, nrow = 1,
widths = c(1,1)
),
ggarrange(
g_models,
labels = c("c"),
ncol = 1, nrow = 1,
widths = 1
),
nrow = 2, ncol = 1,
heights = c(1,1)
)
fig_3births = read_feather(path = str_c(IO$out_Rdata, "simulated_births.feather"))
births_STL = read_feather(path = str_c(IO$out_Rdata, "simulated_births_seasonal_trends.feather"))
opt_par_df = read_feather(path = str_c(IO$out_Rdata, "optimal_parameters_and_AIC.feather"))
ca_levels = dict$country_area$country_area
births = births %>% mutate(country_area = factor(country_area, levels = ca_levels))
births_STL = births_STL %>% mutate(country_area = factor(country_area, levels = ca_levels))
opt_par_df = opt_par_df %>% mutate(country_area = factor(country_area, levels = ca_levels))
ca_wrapped_levels = dict$country_area$country_area %>% str_replace(.," - ","\n")
births = births %>%
mutate(country_area_wrapped = str_replace(country_area," - ","\n") %>%
factor(., levels = ca_wrapped_levels))
opt_par_df = opt_par_df %>% group_by(country_area, BC, sex_type) %>%
mutate(AIC_best_model = min(AIC),
dAIC = AIC - AIC_best_model,
country_area_wrapped = str_replace(country_area," - ","\n") %>%
factor(., levels = ca_wrapped_levels))g_ts = ggplot(births %>% filter(sex_type == "unprot_sex", BC == "all"), aes(x = year_month, y = sim_births, col = model))
g_ts = g_ts +
geom_line(aes(y = births), col = "black", size = 0.7)+
geom_line(size = 0.5)+
geom_text(data = opt_par_df %>% filter(sex_type == "unprot_sex", BC == "all"),
aes(x = Inf, y = Inf, col = model,
label = str_c(AIC %>% round(),ifelse(dAIC > 0, str_c("\n+",dAIC %>% round()),""))),
hjust = 1, vjust = 1, size = 3, lineheight = 0.75)+
guides(col = FALSE)+
xlab("Date")+ylab("")+
facet_grid(country_area_wrapped ~ model, scale = "free" ,switch = "y")+
theme(strip.placement = "outside",
strip.text.y.left = element_text(angle = 0, hjust = 0.5),
strip.background.y = element_rect(fill = "gray90", color = NA))+
ggtitle("Measured vs. simulated births")
g_tsg_st = ggplot(births_STL, aes(x = month, y = sim_births_seasonal, col = model))
g_st = g_st +
geom_hline(yintercept = 0, col = "gray80")+
geom_line(aes(y = births_seasonal), col = "black", size = 1)+
geom_line()+
# geom_text(data = opt_par_df %>% filter(BC == "all", sex_type == "unprot_sex"),
# aes(x = Inf, y = Inf, col = model,
# label = str_c(AIC %>% round(),ifelse(dAIC > 0, str_c("\n+",dAIC %>% round()),""))),
# hjust = 1, vjust = 1, size = 3, lineheight = 0.75)+
scale_x_continuous(breaks = seq(0,12,by = 3))+
guides(col = FALSE)+
xlab("Calendar months")+ylab("Seasonal trends")+
facet_grid(country_area ~ model, scale = "free_y")+
theme(strip.text.y = element_blank())+
ggtitle("Seasonal trends")
g_stg_aic = ggplot(opt_par_df %>% filter(BC == "all", sex_type == "unprot_sex"),
aes(x = model, y = country_area, fill = dAIC))
g_aic = g_aic +
geom_tile()+
scale_fill_gradient(name = "", low = hsv(0.58,0.4,1), high = hsv(0.02,0.9,0.9))+
xlab("")+ylab("")+
#guides(fill = FALSE)+
facet_grid(country_area ~ model, scale = "free")+
ggtitle("Diff. with best AIC")+
theme(axis.text = element_blank(),
strip.text.y = element_blank())
g_aicg_best_aic = ggplot(opt_par_df %>% filter(BC == "all", sex_type == "unprot_sex", model == "A"),
aes(x = country_area, y = AIC_best_model))
g_best_aic = g_best_aic +
geom_bar(stat = "identity")+ coord_flip()+
xlab("")+
facet_grid(country_area ~ model, scale = "free")+
theme(axis.text = element_blank(), strip.text.x = element_blank())
g_best_aicopt_par_df = opt_par_df %>%
mutate(country_area_col = country_area_col %>% factor(., levels = dict$country_area$country_area_col),
Amplitude = 100*alpha,
age_cat_this_ca = get_age_cat(country_area = country_area))
DF = opt_par_df %>%
filter(model != "A",
sex_type == "unprot_sex",
BC == "all",
age_cat == age_cat_this_ca)
g_F = ggplot(DF, aes(x = Tp, y = Amplitude, col = country_area_col))
g_F = g_F +
geom_segment(aes(xend = Tp, yend = 0, linetype = model))+
geom_point()+
scale_color_identity(guide = "legend",
labels = ca_wrapped_levels,
name = "location" )+
scale_x_continuous(limits = c(0,1),breaks = seq(0,3/4,by = 1/4), labels = c("Jan","Apr","Jul","Oct"))+
ylab("Amplitude (% from mean)")+ xlab("Fertility peak time")+
coord_polar(theta = "x", start = 0)+
ggtitle("Calendar year")
g_FDF = DF %>%
mutate(
country = str_split_fixed(country_area, " - ",2)[,1],
Tp_seas = (Tp-10/365 - 0.5*(country == "Brazil"))%%1)
g_F_seas = ggplot(DF, aes(x = Tp_seas, y = Amplitude, col = country_area_col))
g_F_seas = g_F_seas +
geom_segment(aes(xend = Tp_seas, yend = 0, linetype = model))+
geom_point()+
scale_color_identity(guide = "legend",
labels = ca_wrapped_levels,
name = "location" )+
scale_y_continuous(position = "right")+
scale_x_continuous(limits = c(0,1),breaks = seq(0,3/4,by = 1/4), labels = c("Winter\nSolstice","Spring\nEquinox","Summer\nSolstice","Autumn\nEquinox"))+
ylab("Amplitude (% from mean)")+ xlab("Fertility peak time")+
coord_polar(theta = "x", start = 0, clip = "off")+
ggtitle("Solar year")
g_F_seasdf_sine_curves = expand.grid(t = seq(0,12,by = 0.1), country_area = DF$country_area)
df_sine_curves = left_join(df_sine_curves,
DF %>%
filter(sex_type == "unprot_sex", BC == "all", model == "C") %>%
select(country_area, country_area_col, country_area_wrapped, alpha, Amplitude, Tp, Tp_seas) %>%
ungroup(),
by = "country_area")## Adding missing grouping variables: `BC`, `sex_type`
df_sine_curves = df_sine_curves %>%
mutate(
Fertility_calendar = 100 + Amplitude * cos(t/12*2*pi - Tp*2*pi),
Fertility_solar = 100 + Amplitude * cos(t/12*2*pi - Tp_seas*2*pi)
)
g_fert = ggplot(df_sine_curves, aes(x = t, y = Fertility_calendar, col = country_area_col))+
geom_hline(yintercept = 100, col = "gray75")+
geom_line(size = 0.8)+
scale_color_identity()+
scale_x_continuous(breaks = seq(0,12,by = 3), minor_breaks = 0:12, labels = c("Jan","Apr","Jul","Oct","Jan"))+
xlab("Calendar months")+
ylab("Relative changes in fertility (%)")
g_fertg_fert_solar = ggplot(df_sine_curves, aes(x = t, y = Fertility_solar, col = country_area_col))+
geom_hline(yintercept = 100, col = "gray75")+
geom_line(size = 0.8)+
scale_color_identity()+
scale_x_continuous(breaks = seq(0,12,by = 3), minor_breaks = 0:12,
labels = c("Winter\nSolstice","Spring\nEquinox","Summer\nSolstice","Autumn\nEquinox","Winter\nSolstice"))+
xlab("Solar time")+
ylab("Relative changes in fertility (%)")
g_fert_solarfig_4 = ggarrange(
ggarrange(
g_ts,
g_st,
labels = c("a","b","c"),
ncol = 2, nrow = 1,
widths = c(3,1)
),
ggarrange(
g_F + theme(plot.title = element_text(hjust = 0.5), legend.position = "none"),
g_F_seas + theme(plot.title = element_text(hjust = 0.5), axis.title.y = element_blank(), axis.text.y = element_blank()),
ggplot(),
labels = c("e","f",""),
ncol = 3, nrow = 1,
widths = c(1.2,1.8,2), heights = 1
),
nrow = 2, ncol = 1,
heights = c(2,1)
)
fig_4fig_4_plotlist = list()
ordered_ca = c("Brazil - Central-West","United States - California","Brazil - Northeast","United States - Northeast","France","United Kingdom")
model_colors = c("turquoise3","sienna3","maroon3")
# for each country, we put together the time-series and the seasonal trend
for(ca in ordered_ca){
cat(ca,"\n")
age_cat_this_ca = get_age_cat(country_area = ca)
# time-series
this_ca_births = births %>%
filter(country_area == ca,
age_cat == age_cat_this_ca,
BC == "all",
sex_type == "unprot_sex")
this_ca_opt_par_fd = opt_par_df %>%
filter(country_area == ca, age_cat == age_cat_this_ca, sex_type == "unprot_sex", BC == "all") %>%
mutate(dAIC2 = AIC - min(AIC))
g_ts_ca = ggplot(this_ca_births, aes(x = year_month))
g_ts_ca = g_ts_ca +
geom_line(aes(y = births/1000), size = 0.7)+
geom_line(aes(y = sim_births/1000, col = model))+
geom_text(data = this_ca_opt_par_fd,
aes(x = -Inf, y = Inf, col = model,
label = str_c(
model, " : ",
AIC %>% round()#,
#ifelse(dAIC2 > 0, str_c("\n+",dAIC2 %>% round()),"")
)
),
hjust = 0, vjust = 1, size = 4, lineheight = 0.75, fontface = 2)+
scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
scale_x_continuous(breaks = seq(0,2020, by = 5), minor_breaks = 0:2020, expand = expansion(mult = c(0.05, 0))) +
scale_color_manual(values = model_colors)+
facet_grid(model ~ ., labeller = label_both)+
theme(strip.text.y = element_blank())+
guides(col = FALSE)+
xlab("Date")+
ylab("Births [thousands]")+
ggtitle(ca)
# seasonal trend
this_cat_births_STL = births_STL %>%
filter(country_area == ca,
BC == "all",
sex_type == "unprot_sex")
g_stl_ca = ggplot(this_cat_births_STL, aes(x = month))
g_stl_ca = g_stl_ca +
geom_hline(yintercept = 0, col = "gray80")+
geom_line(aes(y = births_seasonal/1000), col = "black", size = 1)+
geom_line(aes(y = sim_births_seasonal/1000, col = model))+
scale_x_continuous(breaks = seq(0,12,by = 3))+
scale_color_manual(values = model_colors)+
guides(col = FALSE)+
xlab("Calendar months")+ylab("Seasonal trends [thousands]")+
facet_grid(model ~ ., labeller = label_both)+
theme(strip.text.y = element_blank())
g_ca = cowplot::plot_grid(g_ts_ca, g_stl_ca, nrow = 1, rel_widths = c(3,1), align = "h")
fig_4_plotlist[[ca]] = g_ca
}## Brazil - Central-West
## United States - California
## Brazil - Northeast
## United States - Northeast
## France
## United Kingdom
fig_4_time_series = cowplot::plot_grid(plotlist = fig_4_plotlist, nrow = 3, ncol = 2, labels = letters[1:6])fig_4_v2 = ggarrange(
fig_4_time_series,
ggarrange(
g_F + theme(plot.title = element_text(hjust = 0.5), legend.position = "none"),
g_F_seas + theme(plot.title = element_text(hjust = 0.5), axis.title.y = element_blank(), axis.text.y = element_blank(),
legend.key.height = unit(22,"pt"), legend.spacing.y = unit(0,"pt")),
g_fert_curves,
labels = c("g","h","i"),
ncol = 3, nrow = 1,
widths = c(1.2,1.8,3), heights = 1
),
nrow = 2, ncol = 1,
heights = c(3,1)
)
fig_4_v2Figures_path = "../Figures Tables Media/Figures/"
scale = 2
ggsave(plot = fig_1, filename = str_c(Figures_path, "F1.pdf"),
width = 17.5, height = 12, units = "cm", scale = scale)
ggsave(plot = fig_2, filename = str_c(Figures_path, "F2.pdf"),
width = 17.5, height = 10, units = "cm", scale = scale)
ggsave(plot = fig_3, filename = str_c(Figures_path, "F3.pdf"),
width = 8, height = 7, units = "cm", scale = scale)
#ggsave(plot = fig_4, filename = str_c(Figures_path, "F4.pdf"),
# width = 17.5, height = 13, units = "cm", scale = scale)
ggsave(plot = fig_4_v2, filename = str_c(Figures_path, "F4.pdf"),
width = 17.5, height = 17, units = "cm", scale = scale)